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1. INTRODUCTION 


A method of orbit determination is investigated which 
employs Picard iteration and Chebyshev series. The method 
is applied to the problem of determining the orbit of an 
earth satellite from range and range-rate observations 
contaminated by noise. The method is shown to be readily 
applicable and to possess linear convergence. 



2. MATHEMATICAL ANALYSIS OF THE PROBLEM 


2.1 FORMULATION OF THE PROBLEM 

Consider the mathematical model of a dynamical system 
subject to observation given by 




X = f(x, 

t) 

(1) 

II 

t) 

(2) 


where 

X denotes an n- dimensional state vector; 

f, an n^vector function of the state; 
y, an ra-vector of observations; 

g, an m- vector function of state variables; and 
t, the time. 

The dot denotes differentiation with respect to time. 

Let the observations be made at discrete time points 
t^ (i=0,l,2,. . . N) in the interval (t^, t^^) . Let these 

observations be denoted by 

y(t^), 1=0, 1,2,. . . N (3) 


In the above model, the functions f and g are nonlinear 
in general and the observations are contaminated by measure- 
ment errors. The problem of state estimation can now be 
stated as follows: 

Find the solution x(t| of the differential system (1) 
such that 


U = 


g (x (t^) , 


t.) - 
1 


y(t,) 


N 

2 

i=0 


2 


(4) 



is a minimum. This, in effect, is to find from the possible 
solutions of the system (1), the solution that satisfies 
the least square criterion with the given observations. If 
the system has a unique solution for a given set of initial 
conditions, then the problem is to determine the set of 
initial conditions, c, for which the solution minimizes the 
function U. 

2.2 PRESENT METHOD 0E SOLUTION 

The present work comprises a study of Picard iteration [1] 
as a method for the solution of the estimation problem. Picard 
iteration has been used to solve the initial value problem [2] 
and the two point boundary value problem [3, 4]. In this 
case the integration constants are determined such that the 
solution satisfies the boundary conditions at each iteration. 
The convergence criteria have been established for this pro^t 
cess. The present work is a natural extension of this method 
to the estimation problem. The technique and the method of 
implementation are discussed in the following sections. 

2.3 PICARD ITERATION 

Picard iteration successively approximates the solution 
of the differential equation (1) using the following 

x^ = f t) (5) 

where i denotes the iteration number and where the initial 
(guessed) solution or zeroth approximation is denoted by 
x^(t). This procedure reduces the solution of the differ- 
ential equation to a sequence of simple integrations. It is 
also seen that at the end of each iteration the integration 
constants c are to be determined. In the case of boundary 
value problems, the constants are chosen such that the 
given boundary conditions are satisfied by the new approxi- 
mate, solution. 



For the present problem, the constants are determined 

so that the solution (t) minimizes the function U given 

by Equation 4. With these new values for the constants and 

the corresponding new approximate solution, the iteration 

process continues to determine the next approximation x. (t), 
^ * . i+1 ' 

It is expected that the iteration procedure converges 

linearly because of its similarity to the classical Picard 
iteration. 

2.4 MINIMIZATION OF THE RESIDUALS 


The minimization of the function U (or the sum of the 
squares of the residuals) can be carried out using any of 
the conventional schemes to solve the function minimization 
problem. For instance, the method of steepest descent could 
be used for this purpose [5]. On the other hand, at the 
extremum, the function satisfies the following 


^c 


N 

2 

i=0 


^x^ Sq 


Sx 


k ) 





( 6 ) 


where 



n X 


1 matrix 


dxj 

^c 


is an n X n matrix 


This is a necessary condition for the extremum of the function 
and gives n algebraic equations for the n unknowns in c. 

These equations are in general nonlinear in c. However, if 
the observations are linear functions of the state, it is 
evident from Equation (6) that this equation set is linear 
in c (as x is a linear function of c). 



The problem ot determining c now reduces to the problem 
of solving the set of nonlinear algebraic equations (6). 

These equations can be solved using Newton method for 
finding the zeros of a function [5, 6J. The derivative matrix 
for L = VU can be written as follows; 


dL' 

be 


= 2 


Sx?' 

k 


X 

/O 

\ 


^c 

axj^ 


J 


+ (g - y(t,)) — 

^ Sc 


Sc Sx, 


This matrix can be used to compute the new value of c using 
the following: 

-1 


c , = c . TV 
J+1 J 


5 ^ 


L . where 
J 


( 8 ) 



is an n X n matrix and 


L. 

J 


corresponds toe.. 


2.5 POLYNOMIAL REPRESENTATION 

In order to implement Picard iteration, a method of 
evaluating the integral (5) must be chosen. The initial 
guessed solution, x (t), and the subsequent iterates, x, (t), 
can be approximated by a polynomial of k-th degree. A poly- 
nomial representation, although other forms are possible, 
has the advantage that it can be integrated easily. 

For a given degree k, the polynomial P, (t) of best 
approximation to the function 0(t) defined on the interval 
(“1, 1) minimizes the norm 

max [4>(r) - P^(t)] (9) 

Polynomials which minimize this norm are said to satisfy the 
rainimax principle [7]. For many functions, a truncated 



series of Chebyshev polynomials is very close to the poly- 
nomial of best approximation [7]. 

For this study , a se;ries/ of Chebyshev -polynomials, is used to 
represent the function Xj^(t). The Chebyshev polynomial of 
k-th degree is given by - 

Tj^(t) ” cos (k arccos (r)) ; “1 ^ t ^ +1 (10) 


These polynomials can be generated from the recurrence 
relations: 




T (t) = 1; T, (t) = T 
O 1 


( 11 ) 


Picard iteration can easily be implemented using the 
orthogonality properties given in Refer;e,nce ' 12. The inter- 
val (t^, t^) is mapped onto the uniform interval (1, -1) by 
the transformation 

T = 1 - 2t/(tj - t^) (12) 

To take advantage of the orthogonality properties, the func^ 
tions are evaluated at the special points given by 


= cos (i 7 t/K) (13) 

where K is degree of polynomial representation. 

The forcing function f(x^s t) as given in Equation 5 
is represented by 

K 

l(x,, t) = 2b. T. (r) 

J=0 J J 


( 14 ) 



The constants b* can be evaluated by knowing the function 
values of f at the special points t j • The integration of 
these Equations (14) can easily be done using the integration 
relations [7J obtaining the polynomial representation for 

Xi^l(T). 



3. APPLICATION TO ORBIT DETERMINATION PROBLEM 


The technique discussed above is employed to solve 
the problem of orbit determination with range and range-rate 
observations from earth-bound tracking stations. The model 
considered and the results are discussed in the following 
sections. All computations were performed on an IBM 360/^65. 

3.1 EQUATIONS OP MOTION 

The equations of motion are formulated in an inertial 
reference frame. This reference :f;rame:fhas: its; or^ 
center of the earth. The X-Y plane coincides with the equa- 
torial plane of the earth with the X-axis pointing to the 
first point of Aries. The Z-axis is in the direction of the 
north pole and the Y-axis forms a right-handed triad with 
the other two axes (see Figure 1) . 

The equations of motion of a satellite, considered 
to be a mass point moving in the central gravity field of 
the earth, in the inertial system of co-ordinates are given 
by [8J 

X = -|J.X/R^ 

.. O 

Y = -p,Y/R'^ 

Z = -M-Z/R^ 

where 
R = (X^ + 

li, = gravitational constant (3.986 x 10^ km^/sec^) 


(15) 

(16) : 
(17) 





3.2 OBSERVATION MODEL 


The earth is modeled as a sphere rotating with an 
angular velocity of 360 deg/day about the Z-axis. The loca- 
tion of the tracking station is specified by giving the 
latitude and the longitude of the station measured with 
respect to the inertial co-ordinate system (XYZ) as shown 
in Figure 1 at the instant t^. 

Two types, of «^obser vat ions are made by each station! 

(1) the range (the distance between the tracking station and 
the satellite) and (2) the range-rate (the time rate of the 


range at the time of observation) 10]. In a practical 

case the range and range-rate are obtained by the radar 


network. These observations are denoted by 
' ^ (j) ' 

A J = 1, 2, . . . S 

= ^(0) i = 0,l,2,...N 

J 


(18) 


Here, the superscript j refers to the station number with 
the total number of stations considered being S. 

The functional form of the observation vector g is 
yet to be given. Let (t), (t) and Z ^ (t) be the 

co-ordinates of the j-th tracking station at time t. Then 
the observation vector as a function of the state variables 
can immediately be written as 


g(t) = 


R^'^^ (t) 

R^^^(t) 


where 

={[X(t) - +[Y(t) - 

+ [Z(t) - 


(19) 


( 20 ) 



( 21 ) 


=-^{[X(t) - [X(t) - X^j^(t)] 

R 

. + [Y(t) - [Y(t) - y^j^(t)] 

+ [Z(t) - [Z(t) - Z^j^(t)]} 

From the above it is seen that both the observations 
and the equations of motion are nonlinear. 


3.3 DETAILS OF NUMERICAL SOLUTION 


The Equations 15 to 17 representing the system dynamics 
are second-border nonlinear differential equations. These 
can be written as a set of six first order nonlinear equations, 
one for each of the state variables (viz, x, Y, Z, X, Y, Z) 
in the form given by Equation 1. The state variables are 
represented by a series of Chebyshev polynomials with 16 
terms, the number of terms being chosen based on the recom- 
mendations of the earlier studies [1, 11] on the representation 
of the orbits using Chebyshev polynomials. The computation 
of the forces is performed only at the selected points as 
given by Equations 12 and 13 in order to maintain the ortho- 
gonality property of the polynomials. The integration of 
the second-order system gives rise to two vector constants 
of integration (each with three elements). Let these six 
constants be represented by c. The function to be minimized 
for the problem is given by 


S 

U = Z 

j=l 


2 ) 
i=0 


- R, 




( 22 ) 


The function given by Equation 22 is minimized with respect 
to c by finding the solution of the system of equations 
given by 



0 


(23) 



Newton's method is used to find the zeros of VU. No major 
difficulties are encountered in implementing this scheme. 

It is observed that U decreases as the solution is 
approached. 

3.4 NUMERICAL RESULTS AND DISCUSSIONS 

To demonstrate the validity of the method and the 
scheme of implementation two examples are considered. The 
first one is the determination of the state of a satellite 
moving in a circular orbit. The orbital elements are given 
in Table I. The observations are made in the interval 
(0, 1513.1) seconds. This corresponds to one- fourth of a 
revolution in the orbit. The tracking station data (viz, 
the location and the number of observations) are also pre- 
sented in Table I. It should be noted that the observations 
are not simultaneous except for the first. This does not in 
any way constrain the applicability of this method for a 
general set of observations. The observations for both 

examples presented here are simulated on the computer using 

th 

a random number generator. At the i iteration the error is 
E = [(X^ - X)^ - (Y. - Y)^ + (Z^ - (24) 

The initial guessed solution for this problem is in 
error by about 72 km throughout the interval considered. In 
Figure 2, the error as a function of time is plotted for the 
first two iterations. From this figure it can be seen that 
although the initial guess is far from the solution, within 
two iterations the error is reduced to less than 5 km. To 
have a clear pictorial representation of the convergence of 
the process, the logarithm of the absolute error is pre- 
sented in Figure 3, as a function of the iteration number. 
The error at selected time points is plotted in order to 



TABLE I 

STATION AND ORBIT DATA (EXAMPLE I) 





Station Number 


.. . 

parameters 

1 

2 

3 


Latitude (deg) 

Longitude (deg) 

Number of observations 

Interval between observations 
(sec) 

Serai-major axis (km) 
Eccentricity 
Inclination (deg) 

Longitude of ascending node 
Argument of perigee 


18.0 


12.0 

ID.O 

0.0 


28.0 

14.0 

10 


20 

30 

168.0 


79.0 

52.0 


a 

= 7178. 145 



e 

= 0.0 



i 

= 20 



= 0 . 0 
CO = 0.0 



Error (km) 



Figure 2, Error as a function^-of time for. different 
iterations (circular orbit) c 


Log of AbS' 



depict the convergence of the overrrall solution. It is 
apparent from this that the error redaction is linear as the 
solution approaches the true solution. 

The convergence criterion is that the error be less 
than a meter. The iteration is terminated when this cri- 
terion is satisfied. «?For this problem, the convergence of 
the process is attained at the end of nine iterations. The 
time for the computation is 21.6 seconds. 

The second problem considered is the case of a 
satellite moving in an elliptic orbit. The orbit has an 
eccentricity of 0.0557 with a perigee height of 400 km. 

Here, the arc considered for the solution is from t^ = 0 
(perigee) to ™ 1513.1 seconds. The circular orbit of the 
first example and this ecdentric orbit have the same semi- 
major axes, and therefore, the same orbital periods. This 
example allows a study of the effect, if any, of eccentricity 
on the method of solution. The tracking stations and the 
times of observation are unaltered from those of the first 
problem. The tracking station data and the orbit data are 
presented in Table II. 

Figure 4 contains the graphical representation of the 
errors in the solution for the first two iterations. The 
error of about 70 km for the guessed solution is reduced 
to about 10 km in two iterations. The logarithm of the 
absolute error is plotted in Figure 5 as a function of the 
iteration number. Here again, it is evident that the con- 
vergence is linear. The process converges, for this problem, 
in 11 iterations to a solution which is in error by about a 
meter. 

In summary, it is jseen that the method is applicable 
to the problem of orbit determination and that the conver- 
gence is linear as expected. The comparison of the two 
problems considered here is presented in Table III. It is 
seen that the elliptic orbit takes more iterations. 



TABLE II 


STATION AND ORBIT DATA (EXAMPLE II) 


Parameters 

'• 1 

Station Number 
2 

3 

Latitude (deg) 

18.0 

12.0 

10.0 

Longitude (deg) 

0.0 

28.0 

14.0 

Number; ^of observations 

10 

20 

30 

‘ t 

Interval between observations 




(sec) 

168.0 

79. 0 

52.0 

Semi-major axis (km) 


a = 7178. 145 


Eccentricity 


e = 0.0557 


Inclination (deg) 


i = 20 


Longitude of ascending node 


= 0. 0 


Argument of perigee 


00 = 0.0 





Figure 4 . Error as a function of time for different 
iterations (elliptic orbit) o 


Log of Absolute Error 




TABLE III 

COMPARISON OF THE TWO EXAMPLES 



Parameter 

Example I 
(Circular Orbit) 

Example II 
(Elliptic Orbit) 

Arc length in time (sec) 

1513. 1 

1513. 1 

Number of observations 

50 

50 

Eccentricity 

0.0 

0.0557 

Initial error in state 
variables 

1.0% 

1.0% 

Number of iterations for 
convergence 

9 

11 

Computer time (sec) 

21.6 

28. 9 


o 




4. CONCLUSIONS 


An iterative scheme using Picard iteration has been 
investigated for the solution of the state estimation pro- 
blem with discrete observations. The applicability of this 
scheme to practical problems has been demonstrated [12] 
using as an example the problem of orbit determination of 
an earth satellite with range and range-rate observations 
from earth-bound tracking stations. 

Unlike some of the more commonly usedi^methods , this 
method does not require the formulation or the solution of 
the linear perturbation equations. Prom the examples con- 
sidered, it is seen the method has linear convergence. 

This scheme can be extended easily to a general problem 
or orbit determination (including the effects of drag, oblate- 
ness of the earth, gravity fields of other bodies, etc.). 
However, further work is necessary to study the convergence 
of Chebyshev series, and the number of terms needed to 
represent the solution over a given interval. 

A sequential algorithm for estimating the orbit using 
Chebyshev series is presently being investigated for use in 
a real time orbit determination environment. 
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